Phase diagram of H2 adsorbed on graphene 
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The phase diagram of the first layer of H2 adsorbed on top of a single graphene sheet has been 
calculated by means of a series of diffusion Monte Carlo (DMC) simulations. We have found that, 
as in the case of He, the ground state of molecular hydrogen is a \/3 x \/3 commensurate structure, 
followed, upon a pressure increase, by an incommensurate triangular solid. A striped phase of 
intermediate density was also considered, and found lying on top of the equilibrium curve separating 
both commensurate and incommensurate solids. 
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I. INTRODUCTION 

Graphene is a novel form of carbon, in which the atoms 
are located in the nodes of a honeycomb lattice that 
extends periodically in tvifo dimensions forming a sin- 
gle layerjiS, In this respect, graphene is different from 
graphite, a substrate formed by the superposition of these 
layers to build a complete three dimensional structure. 
Graphene can be obtained as a free standing structure?^ 
or as a single adsorbed layer on top of another substrate 4 
One of the main interests of the scientific community on 
graphene is related to its unusual electron transport prop- 
erties that are determined by the Dirac equation, '^^^ but 
other characteristics of this compound as for instance its 
behavior as adsorber are starting to be payed attention 
to. A recent diffusion Monte Carlo (DMC) calculation 
of *He on graphene indicates that its behavior is quite 
similar to the one on graphite, the main difference being 
the binding energy, lower in the case of a single carbon 
sheet.® It was also found® that the ground state of helium 
on graphene is a ^/S x -\/3 commensurate solid, in agree- 
ment with experimental data on graphitei^ii^ The aim 
of the present work is to calculate the phase diagram of 
H2 adsorbed on graphene at zero temperature using for 
the first time the DMC method. The comparison with 
the case of graphite will show us if there is any signif- 
icant difference between them or graphene is simply a 
weaker binding version of the graphite phase diagram, as 
happens with ^He. 



II. METHOD 

Our study is based on the DMC method, basically be- 
cause this technique allow us to obtain the correct ground 
state for a given system of bosons. ^^ This is the case here, 
since we will consider only para-H2, the ground state of 
the hydrogen molecule. An important ingredient of any 
DMC calculation is the trial wave function used for im- 
portance sampling. This function collects basic informa- 
tion about the system that is known a priori, and can 



be considered a reasonable approximation (in the varia- 
tional sense) to its ground state. In this work, we use as 
a intial trial wave function 
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that depends on the coordinates of the hydrogen 
molecules ri, r2, ..., tn and the position of the Car- 
bon atoms vj of the substrate. The first term in Eq. 
([1]) is a two-body Jastrow function depending on the 
H2 intermolecular distances r^j, with optimal parameter 
^Ha-Ha = 3.195 A. The second term is another Jastrow 
function that takes into account all the individual C-H2 
interactions, the optimal value of 6c-H2 being 2.3 A. Fi- 
nally, the third term is a product of one-body Gaussians 
that depend only on the z coordinate of each molecule 
and whose function is to localize the molecules near the 
zq value where the binding energy formed by the sum- 
ming up of all the carbon-hydrogen interaction in the 
graphene-adsorbate system is larger. It depends on two 
parameters, whose optimal values are a = 3.06 A~^ and 
Zq = 2.9 A. All the parameters in the trial wave function 
([T]) were variationally optimized for a liquid phase at den- 
sity 0.0068 A^^ and their slight density dependence was 
neglected. 

In the case of the different solid phases considered in 
this work, the trial function above ([T]) was multiplied by 
another set of Gaussians of the form 



JJexp{-c[(Xi - Xsite)^ + iVi - J/sitc)^]} 
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i.e., each particle was limited to be in a region around its 
corresponding Xsitcl/sito coordinates (Nosanow- Jastrow 
model). These coordinates corresponded to the crystal- 
lographic positions of the particular solid considered, and 
the parameter c was variationally optimized for each lat- 
tice type. For instance, for all the commensurate phases 



TABLE I: Equilibrium density po and energy per H2 molecule at equilibrium eo for the liquid phase. Graphene (2D) and 
graphite (2D) indicate the results after the subtraction of their respective energies in the infinite dilution limit. 2D are DMC 
results for a strictly 2D H2 system with the same intermolecular potential (Ref. la). The spinodal densities ps for the same 
systems are also shown. 



Graphene 



Graphite Graphene (2D) Graphite (2D) 



2D 



e(i[K) -451.88 ±0.03 -503.0 ±0.1 -20.09 ± 0.07 -20.4 ±0.1 -21.43 ±0.02 
po (A"2) 0.05948 ± 0.00005 0.0593 ± 0.0004 0.05948 ± 0.00005 0.0593 ± 0.0004 0.0633 ± 0.0003 
Ps (A-2) 0.0489 ± 0.0001 0.0486 ± 0.0001 0.0489 ± 0.0001 0.0486 ± 0.0001 



considered in this work, c = 0.61 A ^. For the incom- 
mensurate solid, a linear fit between the values obtained 
at densities 0.1 (c = 1.38 A'^) and 0.08 A'^ (c = 0.61 
A~^) was used. Below the latter density, c was kept 
fixed to 0.61 A~^. We used the Silvera and Goldman 
potentiali^ for the H2-H2 interaction, one of the stan- 
dard potentials in Monte Carlo calculations of para-H2. 
It considers point-like H2 molecules, due to the low ex- 
centricity of the ellipsoid of the real molecule, and it has 
been shown that reproduces quite accurately the equa- 
tion of state of bulk solid H2. The corrugation effects 
induced by the substrate were introduced by taking into 
account all the C-H2 individual interactions, that are 
assumed to be of Lennard- Jones type with parameters 
taken from Ref. [ij. All the above sets of parameters 
were used both for graphene and graphite alike, the only 
difference in their respective simulations being the num- 
ber of graphene sheets considered in the calculation of the 
C-H2 interaction potential. Following the results of Ref. 
@, we modeled graphite by a series of 8 parallel graphene 
layers separated by a distance of 3.35 Afrom each other 
and stacked in the A-B-A-B form characteristic of this 
compound. The presence of additional carbon layers did 
not modify the H2 binding energies to the substrate, the 
results for 8 and 9 carbon sheets being equivalent within 
their error bars. 



(open boxes). Even though the liquid phase will be shown 
to be metastable, it is useful to calculate its properties 
from a methodological point of view, for the sake of com- 
parison with other possible phases and previous calcula- 
tions. In both cases, we used the same simulation cell of 
34.43 X 34.08 A^, with a variable number of H2 molecules 
on top of it to match the densities displayed in the figure. 
Periodic boundary conditions were considered for the x 
and y directions. To aid in the comparison, the results for 
graphite were up-shifted by the difference in the binding 
energies in the infinite dilution limit between graphene 
and graphite, given in the last paragraph. Since the er- 
ror bars are of the size of the symbols, we can see that 
both curves are virtually identical. This is confirmed by 
a look to the results displayed in Table I, that come from 
least-squared fitting with cubic polynomials to the DMC 
data displayed in Fig. 1. The equilibrium densities po in 
graphene and graphite are equal within error bars, and 
the difference between the energies per molecule at equi- 
librium Co is ^ 0.3 K after having considered the shifting 
due to the differences in the substrate. These results 
indicate that in the metastable liquid phase, the main 
difference between graphene and graphite as adsorbents 
is a nearly constant shift in the binding energies, as in 



III. RESULTS 

We started our study by comparing the experimentally 
obtained binding energies ej,s in the infinite dilution limit 
to the ones calculated by means of the DMC technique 
described above. This was done for graphite, since there 
is no experimental data for graphene yet. The results 
are Cf, = —482.7 ± 3.5 K (experimental, Ref. [Ij) versus 
—482. 57 ± 0.06 K (this work). The good agreement be- 
tween both values means that the C-II2 potential used 
here is expected to give accurate results for this system. 
The energy obtained for graphene was —431.79 ± 0.06 K, 
i.e., a difference with the previous case of 50.78 ± 0.08 K. 

Our main goal was to calculate the phase diagram of 
the first layer of H2 on top of both graphene and graphite 
until the experimental density for promotion to the sec- 
ond layer, 0.0937 A~^— and compare the results ob- 
tained. The first step towards this end is reported in 
Fig. 1, where we show the energies per particle of a liq- 
uid phase on top of graphene (black boxes) and graphite 
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FIG. 1: Energy per H2 molecule in the metastable liquid 
phase for graphene (full boxes) and graphite (open boxes). To 
afford an easy comparison, this latter results are up-shifted 
by the difference between the binding energies of a single hy- 
drogen molecule on top of graphene and graphite. 



TABLE II: Density and energy per H2 molecule for the incommensurate triangular solids obtained by least squared fits to the 
simulation results. Graphene (2D) and graphite (2D) indicate the results after the subtraction of their respective energies in 
the infinite dilution limit, as in the previous table. The 2D entries are results for different calculations of pure 2D II2 systems. 



Graphene 



Graphite Graphene (2D) Graphite (2D) 2D (Ref. 16) 2D (Ref. 27) 2D (Ref. 28} 



eoiK) 
Ps (A-^) 0.0606 ± 0.0001 0.0600 ± 0.0001 0.0606 ± 0.0001 0.0606 ± O.OOOf 0.0584 ± 0.0001 



-454.1 ±0.3 



-505.2 ±0.2 



-22.3 ±0.3 



-22.6 ±0.2 



0.0689 ± 0.0005 0.0689 ± 0.0006 0.0689 ± 0.0005 0.0689 ± 0.0006 



-23.453 ± 0.003 
0.0673 ± 0.0002 



-22.1 ±0.1 
0.064 ± 0.01 



-23.25 ±0.05 

0.0668 ± 0.0005 

0.059 ±0.001 



the case of ""^Hcj^ The spinodal densities ps, defined as 
the points in which the derivative of the pressure with 
respect to the surface density equals zero, are also indis- 
tinguishable in both cases. We can compare the results 
obtained for the liquid phase with the ones for a purely 
two dimensional (2D) system. In Ref. [I6|, this was made 
in the same conditions that the ones considered here, i.e., 
with the same intermolecular potential and for T = K, 
as corresponds to a DMC calculation; their results are 
shown in Table I for comparison. We can infer that the 
introduction of both corrugation and movement in a per- 
pendicular direction to the basal plane have the effect 
of decreasing the binding energy of the liquid at zero 
pressure (21.43 for a flat surface versus around 20 K for 
graphene and graphite). 

Experimental data of H2 on top of graphit o^^'"*^^" — in- 
dicate that the liquid phase is metastable with respect to 
both the incommensurate and commensurate solids, as 
in the *He case. This feature is confirmed by the present 
results shown in Fig. 2, where we display the energies 
of H2 on graphene for its liquid (full boxes), commen- 
surate -\/3 X \/3 structure (open circle, density 0.0636 
A^^), and incommensurate triangular solid phases (open 
boxes). The dotted line is a least-squares fit to DMC data 
for the latter phase with a third degree polynomial. From 
Fig. 2, we conclude that the ground state of H2 on top of 



graphene is a -\/3 x -\/3 commensurate (C) solid, as it is on 
graphite, whose phase diagram is similar and not shown 
for simplicity. The binding energies per H2 molecule for 
this C solid are -461.12 ± 0.01 and -512.97 ± 0.02 K 
for graphene and graphite, respectively. Their difference 
(51.85 K) is slightly larger than the one for the binding 
energy in the infinite dilution limit (50.78 K), indicating 
that the additional stabilization due to the extra graphite 
layers is more important for the C solid phase (1.07 ver- 
sus 0.37 K). The energy of the C phase for graphite is 
slightly below the variational results of Novaco?^ for the 
same structure (—510.8 K). Since the binding energy of 
this commensurate phase is always larger than the one 
corresponding to the liquid phase, we conclude that for 
densities smaller than 0.0636 A~^ the system will break 
in H2 patches separated by empty space to produce the 
average density we could be interested in. 

In order to establish the minimum density for which 
the incommensurate (IC) solid is stable we should make 
a double-tangent Maxwell construction between the C 
and IC phases. Since the C phase is defined by a sin- 
gle density, the construction was made by drawing the 
tangent line to the IC equation of state that intersects 
the C point. The result is depicted in Fig. 3 for the 
case of graphene. There, we can see that the limiting 
density for the IC solid in equilibrium with the v3 x v3 
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FIG. 2: Adsorption energy of H2 on graphene. Full boxes, 
liquid metastable phase; open boxes, incommensurate trian- 
gular solid; open circle, commensurate y/S x y/S registered 
solid. 
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FIG. 3: Same than in previous figure but including the data 
for the a phase (full circles). The x-axis represents now the 
inverse of the density. The error bars are of the size of the 
symbols and not displayed for symplicity. 



TABLE III: Energy per H2 molecule for different arrange- 
ments found in the phase diagram of D2 on graphite and its 
corresponding values for an incommensurate triangular solid 
at the same density. All data correspond to graphene. 



Phase p (A"^) 


Cbinding (K) 


eic (K) 


S 0.0789 
e 0.0835 
7 0.0814 


-442.22 ± 0.08 
-440.30 ± 0.04 
-441.24 ±0.06 


-450.88 ±0.1 
-446.47 ±0.1 
-448.76 ±0.1 



structure is 0.077 A~'^ (the inverse of a limiting surface 
per molecule of 13 A^), in excellent agreement with the 
available experimental data for graphitei^rJS between 9 
and 20 K (0.077 ± 0.001 A'^ ^ 1.22 times the density of 
the C phase). The binding energy corresponding to this 
structure is —452.08 K. The lower density for graphite is 
exactly the same, with a binding energy of —503.15 K. 
Obviously, to draw the limits of this transition is only 
possible if we consider corrugation in the graphene and 
graphite structures. The limits for this C-IC transition 
were not calculated quantitatively in any of the previ- 
ous works of H2 on graphite 1^"— There have been also 
other calculations on the equation of state for a trian- 
gular solid in a purely 2D environment, both at zeroi^ 
and finite temperature i^i2^ Their main results are sum- 
marized in Table II and compared to the present results 
for the IC solid. However, in the context of H2 adsorp- 
tion on graphene and graphite these strictly 2D results 
only serve as a check of the quality of our calculations, 
since the equilibrium density of the triangular 2D solid 
is below the IC density limit determined by the Maxwell 
construction (Fig. 3). 

We have also analyzed the possible existence of the so- 
called a phase, or striped domain phase, in adsorbed H2. 
According to experiments in graphite^i^^ it consists of 
strips of the C phase of variable width separated by nar- 
row walls in which the II2 molecules are closer to each 
other (see the phase diagrams in Ref. KLSl ). This phase is 
also present in D2r^ and appears in the phase diagram 
drawn in Ref l25l The a structure can be defined by a 
rectangular cell that contains more or less H2 molecules 
depending on the width of the C domains. In Fig. 3, 
we display as solid circles the cases for 4,6,8,10 and 12 
molecules per simulation cell. It can be seen that they 
are approximately on top of the Maxwell construction 
line that limits the C-IC transition. This result implies 
that we can consider the a phase as intermediate between 
them, with continuous changes from the C to the a phase, 
and from the a to the IC solid, at least at T = K. This 
continuous change is in agreement with the experimental 



results of Cui et al^ which ruled out a first order change 
between the C and the a arrangement. However, our re- 
sults do not allow us to determine if the a phase is a real 
thermodynamical phase or simply a mixture of the C and 
IC phases separated by striped domains. The reason is 
that the simulation results are below but too close to the 
Maxwell construction line for making a definitive com- 
mitment (see Fig. 3). The data for graphite are similarly 
basicaly on top of the corresponding Maxwell construc- 
tion line, and can be obtained from the graphene ones in 
Fig. 3 by applying a downward shift of 51.46 K. 

To complete the study of H2 on top of graphene, we 
checked the existence of three hypothetical phases that 
appear in the experimental^ phase diagram of D2 on 
top of graphite and that have not been experimentally 
observed in H2. Two of them are commensurate: the e 
phase, which is a 4 x 4 structure {p = 0.0835 A~^), and 
the 5 phase, which corresponds to a 5-^/3 x 5-\/3 arrange- 
ment [p — 0.0789 A~^). Both of them are within the den- 
sity limits corresponding to the IC triangular phase. In 
D2 on graphite, there is also an incommensurate oblique 
phase, which is called the 7 phase, whose range of stabil- 
ity is approximately between 0.077 and 0.083 A^^. We 
calculated the corresponding energies per particle for the 
above mentioned commensurate phases and for a single 
density of the 7 phase (0.0814 A~^), and compared them 
to the energies of an IC triangular phase at the same 
densities. The results are displayed in Table HI. The 
main conclusion is that for H2 all these arrangements are 
metastable with respect to the incommensurate solid. 



IV. CONCLUDING REMARKS 

Summarizing, we have studied the phase diagram of 
H2 on top of graphene using the most powerful micro- 
scopic tool at zero temperature (DMC), accurate poten- 
tials, and incorporating explicit C-H2 interactions (fully 
corrugated model). The phase diagram of the first layer 
on graphene is fully determined for the first time. The 
ground state corresponds to a \/3 x -\/3 commensurate 
solid as happens in ''He4 Graphene and graphite show 
basically the same phase diagram, the main difference 
being the adsorption energy, which is ~ 51 K larger in 
graphite. 
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